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[~ — Abstract 

We propose a method for model reduction on a given frequency range, 
*Vj without the use of input and output filter weights. The method uses a 

nonlinear optimization approach to minimize a frequency limited like 
. cost function. An important contribution in the paper is the derivation of 

the gradient of the proposed cost function. The fact that we have a closed 
I I form expression for the gradient and that considerations have been taken 

to make the gradient computationally efficient to compute enables us to 
T— H efficiently use off-the-shelf optimization software to solve the optimization 

^ problem. 

m 
o 

^ 1 Introduction 

CN Given a linear time-invariant (lti) dynamical model, 

x{t) = Ax{t)+Bu{t), 

• • yit) = Cx{t) + T>u{t) 

where A e M"^", B e M"^"', C e M^^" and D € the model reduction 

problem is to find a reduced order model 

d 

Xr{t) = ArXr{t) + BrU(f), 
yr{t) = CrXr{t) + T)rU{t), 
with Ar e W-xn^^ ^ g ^nr-xm^ g j^pxn. g^j^^j jy g j^pxm ^j^j^ < 

where this reduced order model describes the original model well in some met- 
ric. In this paper we are interested in a reduced order model that describes the 
model well on a given frequency range. This is motivated by situations where 
the given model is valid only for a certain frequency range, for example as in 
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[S] where models coming from aerodynamical and structural mechanics compu- 
tations describing a flexible structure are only valid up to a certain frequency. 

For a review of model reduction approaches, both ordinary and frequency- 
weighted, see e.g. |1] and [5]. Some of the most commonly used frequency- 
weighted methods, according to are [TD], [T] and [B], which all use different 
balanced truncation approaches. In many of the frequency-weighted methods 
one has to specify input and output filter weights. In [5] they introduce a method 
which does not need these weighting functions, by introducing frequency-limited 
Gramians. This method can be interpreted as using ideal low-, band- or high- 
pass filters as weights. However, this method has the drawback of not always 
producing stable models. One approach to remedy this has been presented in 
[3], where they introduce a modification of the method in [5], and additionally 
derive an Hao bound for the error. [5] also presents a modification of the method 
from however this method is only applicable to SISO models. Note that these 
methods, whilst having an upper bound on the error, do not minimize an 
explicit measure. 

One of the main contributions in this paper is a method which uses opti- 
mization, and not truncation, to find an ■H2-optimal reduced order model, and, 
where this optimization is only performed over a limited frequency interval. 
Another important contribution is the derivation of the gradient for the cost 
function. The fact that we have a closed form expression for the gradient and 
that considerations have been taken to make the gradient computationally effi- 
cient to compute enables us to efficiently use off-the-shelf optimization software 
to solve the optimization problem. 



2 Frequency Limited Gramians 



The method proposed in this paper uses the idea presented in [2J, in which 
they introduce frequency-limited Gramians. Before introducing the frequency- 
limited Gramians we define the standard Gramians, see in time and fre- 
quency domain, for reference. For a system G, that is stable and described 

by 



G 



x{t) 

y{t) 



Ax{t) + Bu(t) 
Cx{t) + T)u{t) 



(1) 



denoted as G 
defined as 



A 


B 


C 


D 



the observability and controllability Gramians are 



1 

2^ 



poo 1 pOC 

Q= / e^^^C'^Ce'^^dT = — / H*{iy)C'^Cil{iy)di^, 
Jo 27r 



(2a) 



(2b) 



where H(a;) = {iujl — A)^^ and H*(a;) denotes the conjugate transpose of 
H(w). The controllability and observability Gramians satisfy, respectively, the 
Lyapunov equations 

AP + PA"^ + BB"^ = 0, (3a) 
A'rq -f QA + C"^C = 0. (3b) 
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Now we narrow the frequency band, from (—00, 00) to (— w, lu) where uj < 00. 
We define the frequency-Hmited Gramians, see [5], as 



^u = Y I H(i.)BBTH*(i.)di., (4a) 
Q^ = ^f H*(i.)CTCH(i.)di/. (4b) 
These Gramians can be shown to satisfy the following Lyapunov equations, see 

AP^ + P„AT + S^BBT + BBTs* =0, (5a) 

ATq^ + Q^A + S^CTc + CTcS^^O, (5b) 

with 

S^= ^ln{{A + iLjI){A-iujl)-^) (6) 

In [2] they continue by creating a balanced system and performing a balanced 
truncation using the newly defined frequency-limited Gramians. A drawback 
with this method is that since the terms S„BB"r+BB"rS* and S* C"^C+C"^CS„ 
are not guaranteed to be positive definite, stability of the reduced order model 
cannot be guaranteed. There exist a modification to this in [4 where they 
propose a remedy to this. 

Remark 1. By using addition/subtraction of two or more different frequency- 
limited Gramians it is possible to focus on one or more arbitrary frequency 
ranges, e.g., you can construct the frequency-limited controllability Gramian, 
Pji, for the interval w G il = [wi, W2] U [w3,ijJ4] as 

APn + PoAT + SoBBT + BB^S^^ = 0, (7) 
with Sq — St^2 ^ + ^uii ~ ^uj3 ■ 

3 Frequency Limited Model Reduction using Op- 
timization 

The H2-norm of G, in 0, can be expressed as 

||G||^^ =.tr / Ce^^BBTe^^^CTdr (8a) 
^ Jo 

1 r°° 

= — tr / G{ii^)G*{iiy)diy (8b) 

= — tr / CH(i/)BB"rH*(i/)C"rdi/ = trCPC"^ (8c) 
27r J-00 

= tr / BTe^^^C"rCe^^Bdr (8d) 
Jo 

1 f°° 

= — tr / B"rH*(i/)C"rCH(i^)Bdz^ = trB"^QB. (8e) 
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In this paper we introduce a new frequency-limited ?^2-like norm that uses the 
frequency-limited Gramians presented in the previous section, and we denote 
the new measure by ||G||^^ ^, with 

^^tv j^^^G{iv)G*{ii.)dv (9) 

= — tr/ (CH(i/)B + D) (B"rH*(i/)C"r + D"r)dt/ (10) 
27r J _^ 



:trCP^C"^ + 2tr 



CS^B + D 



27r 



D' 



(11) 



One thing that differs from the ordinary H2-norm is that, if we do not include 
an infinite interval in f2, i.e., include w = cxd as the end frequency, then the 
system does not need to be strictly proper. This means that we can, in this 
case, have D 7^ 0. 

The method proposed in this paper is a model reduction method that, given 
a model G, finds a reduced order model, G, that is a good approximation on 
a given frequency interval, e.g. [0,w]. The objective is to minimize the error 
between the given model and the sought reduced order model in a frequency- 
limited ■H2-iiorm, using the frequency-limited Gramians. We formulate the op- 
timization problem 



minimize 
G 



G-G 



= minimize 

'H2,i^ G 



\E\ 



where 



\E\\l =— tr 
I ii«2," 27r 



E{ii')E*{iu)dv. 



Assume that the system E is stable and described by 

E 



Ae 


Be 


Ce 





Given G and G, represented as 



E 



G : 


A 


B 


C 


D 


be realized. 


in 


Ae 


Be 




Ce 


De 





,G 



B 



D 











1) 




(c 




D-D 



(12) 



(13) 



(14) 



(15) 



(16) 



This realization of the error system will later prove beneficial when rewriting 
the optimization problem. Throughout the paper we will assume that the given 
model is stable. 



The cost function for the optimization problem (12 1 can be written as 



\E\\l^^^^=ti'CEPE,.Cl + 2ti- 
= trB];Q£;,^BE-f-2tr 



CeSemBe - 

CeSe,uiBe 



D 



D, 



'2tt 

LO 

'2tt 



D 



(17a) 
(17b) 
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where 



Qe,oj-^E ■ 



- ^Em^e'B'e ' 
^*E,LU^~E^E 



0, 
0, 



with 



Sb^c = ^ In {{Ae + iu;1){Ae - iljI)-^) 



(18a) 
(18b) 

(19) 



In this paper we have also derived a simpler expression for the matrix S^^^, 
compared to what is presented in [2]. 



Lemma 1. For a matrix A that is Hurwitz we have that 



2tt 



ln((A + iwI)(A-iwI)"i) 



Re 



In (-A - iul) 



Proof. See|J] 



(20) 



□ 



Now we want to rewrite the cost function (17 1 to a more computationally 
tractable form. This is done by using the realization given in ( 16 1 and by 
partitioning the Gramians Pe,uj smd Qe,u: as 



XT 



p. 



Qe,. 



and Se Lj sls 







(21) 



(22) 



Pw, Qw,Pw, QwjXtj and satisfy, due to (18 1, the Sylvester and Lyapunov 
equations 

(23a) 
(23b) 
(23c) 
(23d) 
(23e) 
(23f) 

ith 



AP^4 


P.;AT 




hBB^s: 


-0, 


AX„4 


X^AT- 


h S^BB^ 


hBE^s: 


-0, 


AP^^ 


- P..AT - 


h S^BB^ 


hEB^s: 


= 0, 


A^Q^ 


+ QmAh 






-0, 


A^Y^ 


+ Y„A- 


- S^C^C - 


- c^cs^ 


= 0, 


A^Q^ 


+ Q<^ah 






= 0, 



= Re 



27r 



ln(-A - iiol) 



S„ = Re 



27r 



In (^-A - iul^ 



(24) 



Note that P^^ and Q^^ satisfy the Lyapunov equations for the frequency- 
limited controllability and observability Gramians for the given model, and P^^ 
and satisfy the Lyapunov equations for the frequency-limited controllability 
and observability Gramians for the sought model. 
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With the partitioning of Pe,uj and Qe.u it is possible to rewrite ( 17) in two 
alternative forms 



2tr 



CS^B+D 



2^ 



( 



CS ,B + D 



'2i 



\E\ 



tr ( CP,„CT - 2CX,„CT + CP,„CT 



2tr 



CS,,B + D: 



CS,,B + D: 



D' D' 



(25a) 



(25b) 



27r V 27r. 
Remark 2. Note that neither the term B^Q^^B nor the term CP^jC^, which 



are included in the cost function (251, depend on the optimization variables, 
A,B,C and D. Hence, these terms can be excluded from the optimization. 
These are the only terms including P^^ and Q^^ which are the most costly to 
compute. 

When optimizing the frequency-limited 7{2-iiorm using the system matrices 
as optimization variables we have the freedom to choose which elements we 
optimize over, i.e., we can introduce structure in A,B, C and D, as long as we 
can find an A that is Hurwitz. Let us introduce the matrices S^, Sg, and 

which hold the structure of the sought matrices, i.e.. 



if A 



is a free variable; 

i 

otherwise. 



(26) 



We will see in the next section, that due to the element-wise differentiation, this 
structure will be inherited in the gradient. 

The parametrization of the sought system using the full system matrices is 
of course redundant, which leads to a non-unique minimum of the cost function 
in the parameter space. This leads to a singular Hessian matrix. However, this 
is taken care of in most quasi-newton solvers to ensure that the minimum is 
reached in a numerically stable way. 



3.1 Gradient of the Cost Function 

An appealing feature of the proposed nonlinear optimization approach, using 



our proposed ?^2-like measure to solve the problem, is that the equations (251 
are differentiable in the system matrices. A, B, C and D. In addition, the closed 
form expression obtained when differentiating the cost function is expressed in 
the given data (A,B,C and D), the optimization variables (A,B,C and D) 



and solutions to the equations in ( 23 1 



To show this we start by differentiating with respect to B, C and D. First we 



note that neither Q^^ , nor Q^^ in equation ( 25a ) depends on B which means 



that the equation is quadratic in B. Analogous observations can be made with 



equation (25b I and the variable C and similarly with D. Hence, the derivative 
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of the cost function with respect B , C and D becomes 



d\\E"^ 



= 2 (Q.B + Y^B - S^CT (d - d) ) Sg, (27a) 



9\\E\\l.,,^ 



2(CP„-CX<,- (d-d)bTs3) 0S<5, (27b) 



d\\E"^ 



dC 

--2 fcS„B + D;^-CS„B-D:^+ fD-Dl 0S 



dB V 27r 2tt \ J 2tt 



D' 

(27c) 



where represents the Hadamard product of matrices, i.e., element- wise mul- 
tiphcation. 

For the more comphcated case of differentiating with respect to A we observe 
that and Y^^ do depend on A, see the equations in (23 1. The calculations 
of this part of the gradient are lengthy and can be found in |B[ 

The complete gradient becomes 



d\\E\\l^,. 
dA 

dB 

d\\E\\l^,. 
dC 

dl) 

where 



=2(y3x+Q„p) 0S^-2W0S^, (28a) 
^2(ci^B + YlBSlC'' (P--d))qS^, (28b) 



2 (CP^ - CX^ - - Dj B^S^j Sc, (28c) 
- 2 (CS^B + Dw - CS^B - Dw + (d - d) S^,, (28d) 



W=Re ^l(-A- iujl,Y^ , (29a) 

V =C^CP - C"rCX - C"^ - B"^ (29b) 

with the function _L(-, •) being the Frechet derivative of the matrix logarithm, 
see |S]. 

Remark 3. Remember that if uj = oo is included as the end frequency we need 
to have D^; = D — D = 0, i.e., T) constant and with D = D. 

Remark 4. The proposed method can, analogous to what is done in ^T^, easily 
be extended to a method for identifying LPV-models over a limited frequency 
domain. 

Remark 5. By supplying the cost function and its gradient in computationally 
efficient forms this method can be used in any off-the-shelf Quasi-Newton solver. 

Remark 6. By using a stable model, e.g., a model from a Hankel reduction, 
as an initial point in the optimization and using a line-search we can limit the 
search to stable models. 



7 



4 Numerical Examples 



In this section three examples are used to illustrate the applicability of the 
method and to compare it with other methods. In the examples we will use 
three different methods; Truncation of Hankel singular values (will be called 
Hankel), the method proposed in [2] (called Gawronski), the method proposed in 
[Ij (called Mod. Gawronski) and the proposed method (called Prop, method). 
The proposed method and the original and modified Gawronski method will 
take the limited frequency range into account, but Hankel will not. We will 
also compare with the Hankel method where we use an input filter to help that 
method to focus on the frequency interval of interest. The Gawronski method is 
a representative method among frequency- weighted methods, see [1], with the 
benefit of not having to design weighting functions, but with the drawback that 
it cannot guarantee that the resulting model is stable. 

The proposed method uses a cost function which is non-convex, which makes 
it important to use a good initial point. For the examples presented here we 
have used the model obtained by the Gawronski method as an initial point for 
the optimization in the proposed method. If the Gawronski method generates 
an unstable model we use the Hankel method instead. 

To evaluate the different models against each other we will compare the error 
model, G — G, for the given frequency interval using the limited frequency 'H2- 



norm, 



G-G 



G-G 



, and the relative limited frequency H2-norm, 
We will also compute the relative Hoo-norm on the given frequency interval, 

\\G-G\\ 

denoted — rrrrr] — , and the eigenvalue with the largest real part. 

Example 1 (Small illustrative example). This example addresses a small model 
with four states. The model is composed of two second order models in series, 
one with a resonance frequency at uj ~ I and the other ai w = 3. We will limit 
the frequency range to lo G [0, 1.7] to try to only capture the first model. 

1 9 

G - G1G2 = ^ Q + 1 s2 + 0.003,s + 9- ^^^^ 

To help the Hankel method we create a low pass Butterworth filter of order 10 
with a cut off frequency of 1.7, see Figure^ The results from the different 
methods can he seen in Figures [1| and [5| and Table [7] As can he seen in the 
result we are successful in finding a good model for the first model with both 
the proposed method and the Gawronski method. All the reduced order models 
are stable. The Hankel method and the modified Gawronski method captures 
the wrong resonance mode (from our perspective) and fails completely in the 
lower frequency region. The Hankel method on the filtered model finds a model 
with the same resonance frequency, but otherwise has a bad correspondence with 
the true model. The proposed method and the Gawronski method return models 
essentially indiscernible. 

Example 2 (Example 1 in [5]). In this example we reuse Example 1 from J^. 
The model, which is a spring-damper model with three masses, has six states 
and we will reduce the model to three states and limit the frequency interval to 
u) G [1.5,3.2]. To help the Hankel method we create a band pass Butterworth 
filter of order 10 with cut off frequencies of 1.5 and 3.2 rad/s, see Figure^ 
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Magnitude plot for the filter, true and filtered true model 



30 



20 



10 



TJ ■ 



-10 
-20 
-30 
-40 
-50 



- G*W 



- W 

G 



1 1 

' i 



;^o-0.4 ^0-0.2 ^(,0 ;^(,0.2 

Frequency [rad/s] 



10" 



10" 



Figure 1: Magnitude plot of the given model, the filtered model and the filter 
used in Example [T] For a more easy read plot, in color, the reader is referred 
to the digital version of this paper. 
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Magnitude plot for the true and the reduced models 
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Figure 2: Magnitude plot of the given and reduced order models in Example [T] 
for w e [0, 1.7]. The dashed black vertical line denotes uj = 1.7. We see in the 
figure that the proposed method and the Gawronski method finds models which 
are good approximations of the true model on the given frequency interval. For 
a more easy read plot, in color, the reader is referred to the digital version of 
this paper. 
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Magnitude plot for the error models 
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Figure 3: Magnitude plots of the error models in Example [T| for uj G [0,1.7]. 
The dashed black vertical line denotes uj — 1.7. We see in the figure that the 
proposed method and the Gawronski method finds accurate approximations to 
the true model on the given frequency interval. For a more easy read plot, in 
color, the reader is referred to the digital version of this paper. 
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Table 1 : Numerical results for Example 







G-G 




\\g-g\\ 
















Hankel 


1.77e4 


00 


l.Ole+00 


l.OOe+00 


-1.59e-03 


Hankel w. filter 


2.93e+00 


1.67e+00 


2.14e+00 


-4.03e-02 


Mod. Gawronski 


1.77e+00 


l.Ole+00 


l.OOe+00 


-1.51e-03 


Gawronski 


9.14e-02 


5.21e-02 


3.35e-02 


-9.88e-02 


Prop, method 


8.51e-02 


4.85e-02 


3.26e-02 


-9.94e-02 



Magnitude plot for the filter, true and filtered true model 
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Figure 4: Magnitude plot of the given model, the filtered model and the filter 
used in Example [2] For a more easy read plot, in color, the reader is referred 
to the digital version of this paper. 



The results from the different methods can be seen in Figures^ and^ and Table 
The proposed method and the Gawronski method are also in this example 
successful in finding low order models that approximate the given model on the 
given frequency range, and all the reduced order models are stable. The proposed 
method and the Gawronski method captures the correct frequency interval with 
good accuracy. The Hankel method on the filtered model captures the correct 
resonance frequency, however, with the worst overall 'H2-^eo,sure. The other 
two methods misses the relevant frequency interval. Only the proposed method, 
the Gawronski method and the Hankel method with the filtered model are able 
to find the right resonance peak, whilst the modified Gawronski and the Hankel 
method are completely off. 

Example 3 (Aircraft example). The model in this example is a model with 
22 states that describes the longitudinal motion of an aircraft, see J^. We will 
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Magnitude plot for the true and the reduced models 




Figure 5: Magnitude plot of the given and reduced order models in Example 
[2] for u! e [1.5,3.2]. The dashed black vertical lines denotes ui — 1.5 and 3.2. 
We see that proposed method and the Gawronski method tries to capture the 
correct resonance frequency with good accuracy, also the Hankel method on the 
filtered model does this but not as good. The other methods misses the relevant 
frequency interval. For a more easy read plot, in color, the reader is referred to 
the digital version of this paper. 
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Magnitude plot for the error models 
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Figure 6: Magnitude plots of the error models in Example |2] for w G [1.5,3.2]. 
The dashed black vertical lines denotes w = 1.5 and 3.2. The proposed method 
and the Gawronski method finds the best approximations on the relevant fre- 
quency interval. For a more easy read plot, in color, the reader is referred to 
the digital version of this paper. 
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Table 2: Numerical results for Example^ 







G-G 




\\g-g\\ 
















Hankel 


9.27e-01 


3.05e+00 


9.59e-01 


-3.54e-03 


Hankel w. filter 


4.58e-01 


1.50e+00 


1.71e+00 


-2.02e-02 


Mod. Gawronski 


3.10e-01 


1.02e+00 


9.89e-01 


-3.78e-03 


Gawronski 


3.30e-02 


1.09e-01 


4.36e-02 


-2.86e-02 


Prop, method 


2.72e-02 


8.94e-02 


5.34e-02 


-2.86e-02 



Magnitude plot for the filter, true and filtered true model 
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Figure 7: Magnitude plot of the given model, the filtered model and the filter 
used in Example [3] For a more easy read plot, in color, the reader is referred 
to the digital version of this paper. 



reduce this model to 15 states and limit the frequency range to uj £ [0,15]. 
To help the Hankel method we create a low pass Butterworth filter of order 
10 with a cut off frequency of 15 rad/s, see Figure^^ The results from the 
different methods can be seen in Figures^ and^ and Table^ In this example 
the Gawronski method results in a model which at a first glance looks very good, 
however the model is unstable. The other methods capture the given model with 
varying accuracy. Since the Gawronski method generates an unstable model in 
this example, the proposed method is initialized with the model from the Hankel 
reduction instead. The proposed method finds the most accurate model on the 
given frequency range. 

In all three of the above examples we observe that the proposed method 
finds an at least as good model as the method proposed in [2j, but with the 
proposed method we can guarantee that the reduced order model is stable and 
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Magnitude plot for the true and the reduced models 




Figure 8: Magnitude plot of the given and reduced order models in Example [3] 
for ui € [0, 15]. The dashed black vertical line denotes uj = 15. The Gawronski 
method looks to have found the best approximation, however the model found 
by the Gawronski method is unstable. The proposed method finds the best 
stable approximation of the true model. For a more easy read plot, in color, the 
reader is referred to the digital version of this paper. 
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Magnitude plot for the error models 
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Figure 9: Magnitude plots of the error models in Example [3] for oj e [0, 15]. The 
dashed black vertical line denotes uj — 15. The Gawronski method looks to have 
found the approximation with the smallest error, however the model found by 
the Gawronski method is unstable. The proposed method finds the model with 
the lowest error which is stable. For a more easy read plot, in color, the reader 
is referred to the digital version of this paper. 
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Table 3: Numerical results for Example 







G-G 




\\^-^\\u... 












I|G||„,,„ 




Hankel 


5.15e4 


-00 


1.39e-01 


1.16e-01 


-4.85e-03 


Hankel w. filter 


5.10e4 


-01 


1.37e+00 


1.17e+00 


-1.24e-01 


Mod. Gawronski 


5.00e+00 


1.35e-01 


1.30e-01 


-3.19e-03 


Gawronski 










4.67e+00 


Prop, method 


1.68e-01 


4.51e-03 


1.51e-02 


-1.77e-01 



also impose structure in the system matrices. In the first two examples it takes 
less than one second and in the third example about 20 seconds for the proposed 
method to find a reduced order model. 

5 Conclusions 

In this paper we have proposed a method, based on nonlinear optimization, 
that uses the frequency- limited Gramians introduced in and the method 
does not have the drawback of finding unstable models. We use these Grami- 
ans to construct a frequency-limited ■H2-iiorm which describes the cost function 
of the optimization problem. We derive a gradient of the proposed cost func- 
tion which enables us to use off-the-shelve optimization software to solve the 
problem efficiently. The derivation of the method also enables us to impose 
structural constraints, e.g., upper triangular A-matrix, in the system matrices. 
The derivation follows closely the technique in |2j and it is easy to extend the 
method to identify LPV- models. We also presented three examples of different 
sizes and characteristics to show the applicability of the method. The method 
uses a nonlinear optimization approach on a nonconvex problem, hence the un- 
avoidable problem of local minima are present. However, we consider this as a 
two step approach where you use your favorite reduction method as an initial 
point to this method which will refine the solution, as we have shown in the 
examples. 
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A Proof of Lemma [T] 

We want to show that 



27r 



Re 



In (-A - iul) 



(31) 



assuming that A is Hurwitz, i.e., A is a real matrix and if A is an eigenvalue 
to A then we have that Re (A) < 0. We can decompose A as A = VDV~^ 
where D is a diagonal matrix with the eigenvalues to A on the diagonal and 
V is the matrix with the eigenvectors to A as its columns. For a real matrix 
we have that the eigenvectors for real eigenvalues are always real and that the 
eigenvectors for a complex-conjugate pair are also complex-conjugate. We can 
write A as 

"Ai ••• 0' 

A2 

A=VDV-' =[viV2---Vn] . 

.0 ••• A„ 

= ViViXi + U2S2A2 H h VnVnK 

























yl. 



(32) 
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and /(A) as 



/(Ai) 
/(A2) 

/(A) =V/(D)V-i = [v^V2 



='!;i'i;i/(Ai) + V2V2f{>^2) H h w„i)„/(A„) 



/(A„). 









\vj] 
















-I. 



(33) 



one term in equation ( 33 1 and 



For a real eigenvalue, A, in A with corresponding eigenvector, v, represented as 

I and 

/(A) =iln ( ^-i^ ) =iln{-X-iuj)-iln{~X + iuj) (34) 



A — iuj ^ 
we can write this term as 

vv'^ {iln (—A — iuj) — iln (—A + iu) ) = vv~^ (^iln (—A — iuj) + iln (—A — iuj)^ 
= 2vv'^ Re In (- A - iw) ] = 2 Re [vv'^i In (- A - iw) ] =2 Re [wi;"^g( A)] . (35) 



If we now instead look at two terms in ( 33 1 corresponding to a complex-conjugate 



pair of eigenvalues, A and A, with corresponding eigenvectors, v and v, then we 
can write this as 



vv~^ f{\) + vv~^ f{X) ~ vv~^ i In (—A — ito) — vv~^ i In (—A + iuj) 

+ vv^i In (—A — iuj) — vv^i In (— A + iuS^ = vv~^ g{\) + vv~^ g{\)+ 
+ vv 



7Tg(A) + wTg(A) = 2Re [vv'^g{X)] + 2Re 
This means that for a matrix, A, which is Hurwitz, we have 



vv^g{X) . (36) 



/(A) = V/(D)V-1 = ViVifiXi) + V2V2f{X2) + ■■■ + VnVnfiXr.) 

= 2Re [viVig{Xi) + W2W2ff(A2) H h w„i)„5(An)] = 2 Re 5(A) 



i.e., we have 



= — In ((A + iwI)(A - iojl)-'^) = Re 



In (-A - iLoI) 



(37) 



(38) 



B Derivation of the Gradient w.r.t. A 

In this appendix we present the differentiation of the cost function with respect 
to A. 

The cost function is 



(b^Q^B + 2BTy^B + B^Q^B 



2tr 



CS<,B + - f CS^^B + D;^ 



(39) 
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and by looking at the equations 

A^Y^ + Y^A - S* C^C - C^CS^ = 0, 
A^Qc. + Q^A + S*^C'^C + C^CS^ = 0, 



(40a) 
(40b) 



we observe that Q^^ and Y^j depend on A which we need to keep in mind when 



differentiating (|39| with respect to A. Hence 



dA 



tr 2BB ' 



da.; 



BB 



OA 



2C 



becomes 



B D 



D' 



(41) 

where and depend on A via the differentiated versions of the equations 
in (|40|, 



dYl 

— H ^ 

dcLij dciij 



rT _ ^^ZqTq 

■ " da,. 



0, 



^A+^Q +Q ^ 

aaij oaij oaij aaij 



(42a) 
0. (42b) 



To be able to substitute and to something more computationally 

tractable we use the following lemma. 

Lemma 2. // M and N satisfy the Sylvester equations 

AM + MB + C = 0, NA + BN + D = 0, 

then trCN = trDM. 

Studying Lemma |2j the two factors in front of and in (41 ) and the 

structure of the Lyapunov/ Sylvester equations in (42 1, A^ • + • A + * = and 
A^ • + • A + = 0, brings us to the conclusion that, to do the substitution, we 
need to solve two additional Lyapunov/ Sylvester equations, namely 



AX + XA^ 
AP + PA^ 



BB' 



0, 



BB' =0. 



(43a) 
(43b) 



Note that P in ( 43b I is the controllability Gramian for the reduced order model. 
Rewriting (41 1 using Lemma [2] ( [42| ) and (43 1 leads to 



\9\\E\\IJ 


=2tr 


dA 





dA'^ 



da 



2tr 



da. 



- - \ dS* / ^ ^ ^ ^ \ 
Y^X + Q^P j + (C^CP - C^CXj 

B (d^ 



■12 



(44) 



What remains is to rewrite the two last terms in ( 44 1 , which includes and 



as; 

dai. 



Recall the definition of S^^, 



S,„ = Re 



In I —A - iuA 



= Re 



In r(A) 



(45) 
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and differentiate with respect to an element in A, i.e., a.y . This yields 



dai. 



= Re 



27r 



L r(A) 



dr{A) 
dan 









)] 


= Re 








2tt \ 



OA 



(46) 



where L(A,E) is the Frechet derivative of the matrix logarithm, see [S], with 

L(A,E)=/' (i(A-I) + I)~^E(i(A-I)+I)"Mt, (47) 
Jo 

r(A) = - A - iuji. (48) 

The function L(A,E) can be efficiently evaluated using the algorithm by 0. 

By substituting (46) into (44 1 and using (47 1 with the fact that we can 
interchange the tr-operator and the integral we obtain 



\9\\E\\l,,J 


= 2tr 


OA 





da. 



2tr 



2tr 



Y^X+Q^P 



tr 



- 2tr 
dciii 













)*] 


9AT 




Re 





1l(ka),v) 



2 ( Y^X + Q^P ) - 2W 



where 



L r(A),V 



W =Re 



V =C^CP C^CX ( D D I 



(49) 

(50) 
(51) 
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